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ABSTRACT 

Upcoming space-based surveys such as Euclid and WFIRST-AFTA plan to measure Baryonic Acous- 
tic Oscillations (BAOs) in order to study dark energy. These surveys will use IR slitless grism spec- 
troscopy to measure redshifts of a large number of galaxies over a significant redshift range. In this 
paper, we use the WFC3 Infrared Spectroscopic Parallel Survey (WISP) to estimate the expected num- 
ber of Ha emitters observable by these future surveys. WISP is an ongoing Hubble Space Telescope 
slitless spectroscopic survey, covering the 0.8 - 1.65/im wavelength range and allowing the detection 
of Ha emitters up to z ~ 1.5 and [OIII] emitters to z ~ 2.3. We derive the Ha-[0III] bivariate line 
luminosity function for WISP galaxies at z ~ 1 using a maximum likelihood estimator that properly 
accounts for uncertainties in line luminosity measurement, and demonstrate how it can be used to 
derive the Ha luminosity function from exclusively fitting [OIII] data. Using the z ~ 2 [OIII] line lu- 
minosity function, and assuming that the relation between Ha and [OIII] luminosity does not change 
significantly over the redshift range, we predict the Ha number counts at z ~ 2 - the upper end of 
the redshift range of interest for the future surveys. For the redshift range 0.7 < z < 2, we expect 
^3000 galaxies/deg 2 for a flux limit of 3 x 10“ 16 ergs s _1 cm -2 (the proposed depth of Euclid galaxy 
redshift survey) and ^20,000 galaxies/deg 2 for a flux limit of ^10“ 16 ergs s -1 cm -2 (the baseline 
depth of WFIRST galaxy redshift survey). 

Subject headings: galaxies: high-redshift, luminosity functions, number counts 


1. INTRODUCTION 

The origin of dark energy, responsible for the acceler- 
ated expansion of the universe first observed by Riess et 
al. (1998) and Perlmutter et al. (1999), is one of the most 
important unsolved problems in cosmology today, and 
significant effort is being devoted to constrain its prop- 
erties. Dark energy affects both the expansion history of 
the Universe as well as the growth of structures. Both 
effects can be observationally constrained through large 
galaxy redshift surveys, which enable the measurement of 
Baryon Acoustic Oscillations (BAOs; thus, constraining 
the cosmic expansion history) and large scale redshift- 
space distortions (thus, constraining the growth history 
of the large scale structure). The combination of these 
two measurements allows the differentiation between an 
unknown energy component and modification of general 
relativity as the cause of observed cosmic acceleration 
(Guzzo et al. 2008; Wang 2008). 

Upcoming space-based missions, ESA’s Euclid (Lau- 
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reijs et al. 2012) and NASA’s WFIRST-AFTA (Green 
et al. 2012; Dressier et al. 2012; Spergel et al. 2015), 
will perform complementary galaxy redshift surveys to 
map the large-scale structure and its evolution over a 
cosmic time covering the last 10 billion years. Both Eu- 
clid and WFIRST-AFTA will use the Ha A6563 line and 
[OIII] AA4959+5007 doublet to select emission line galax- 
ies as tracers of the large scale structure at 0.7 < 2 (in 

Ha) and 2 < 2 : <2.7 (in [OIII]). The performance of the 
planned missions can be quantified by a figure-of-merit, 
which describes their ability to measure the present value 
and time evolution of the dark energy equation of state. 
The figure-of-merit for a dark energy survey depends on 
the number density of tracer galaxies available at each 
redshift. It is therefore critical to have a reliable and 
sufficiently precise knowledge of the expected number of 
Ha and [OIII] galaxies in the survey volumes. 

Euclid and WFIRST-AFTA will both perform IR slit- 
less spectroscopy of emission line galaxies, in a way sim- 
ilar (scaled by many orders of magnitude in area) to the 
IR spectroscopic surveys that are being conducted with 
the Wide Field Camera 3 (WFC3) on-board the Hub- 
ble Space Telescope (HST). The WFC3 Infrared Spec- 
troscopic Parallel Survey (WISP; Atek et al. 2010) is an 
on-going pure-parallel near-infrared grism spectroscopic 
survey using the WFC3 camera. While covering a sub- 
stantially smaller area, WISP is very similar in many re- 
spects to the planned dark-energy surveys, and thus, can 
be used to test number count predictions, redshift mea- 
surement accuracy, target selection function, as well as 
completeness for the planned surveys. Towards this goal, 
Colbert et al. (2013) predicted the number counts of Ha 
- emitting galaxies in the redshift range 0.3 < z < 1.5. 
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In this paper we extend the work by Colbert et al. (2013) 
and estimate Ha number counts out to z ~ 2. 

Ground-based wide-held narrow-band surveys like 
HiZELS (Geach et al. 2010; Sobral et al. 2009, 2012, 
2013) and NEWFIRM Ha Survey (Ly et al. 2011) have 
been able to measure the Ha luminosity function in the 
redshift range of interest (0.7<2:<2). However, while 
having the advantage of high sensitivity to emission lines 
and covering significant areas in the sky, these surveys 
can only map very narrow redshift ranges. Volume den- 
sities of galaxies can thus be strongly affected by the 
presence of large scale structures in the held. Moreover, 
samples selected with narrow band surveys, without an 
extensive spectroscopic followup, can suffer from contam- 
ination by emission lines at different redshifts (e.g., [OKI] 
and [Oil]; Martin et al. 2008; Henry et al. 2012). Finally, 
even narrow-band surveys with multiple hlters tuned to 
identify multiple emission lines at the same redshifts still 
rely on continuum detections and miss the lowest mass 
galaxies, to which WISP is very sensitive. 

WISP’s grism coverage includes Ha for the redshift 
range 0.3 < z < 1.5 and [Olll] for 0.7 < z < 2.3. Since 
Ha is not directly covered by WISP at z > 1.5, we cannot 
measure the Ha luminosity function up to z^2 explicitly. 
However, we can use the [Olll] coverage to estimate the 
Ha luminosity function and number counts. In order to 
do this, we compute the bivariate Ha-[OlIl] line luminos- 
ity function in the redshift where both lines are visible, 
and use the resulting ht at higher redshifts, where only 
[Olll] -emitters are observable in the WISP data. We also 
introduce a modihed Maximum Likelihood Estimator to 
obtain the best-fit model parameters, which accounts for 
measurement uncertainties in the line luminosity (which 
can substantially affect the shape of the bright end of the 
luminosity function). 

The paper is organized as follows: in Section 2, we 
summarize the new WISP data that we use in this work; 
Section 3 discusses the observed Ha-[OlIl] relation; and 
Section 4 describes the parametrization for the bivari- 
ate luminosity function. In Section 5, we discuss the 
Maximum Likelihood Estimator modified to account for 
uncertainties in the line flux, and the fitting procedure, 
which we use to derive the Ha-[OlIl] bivariate luminos- 
ity function at redshift z ~ 1 and the result is discussed 
in Section 6. Further, in Section 7, we demonstrate the 
ability to recover the Ha luminosity function as well as 
number counts from fitting only the [Olll] data at red- 
shift z ~ 1. Lastly, we fit the [Olll] luminosity func- 
tion and use it to derive the Ha luminosity function and 
number counts at redshift z ~ 2 in Section 8, along with 
the final number count estimates for the upcoming dark- 
energy surveys. 

Throughout this paper, we assume standard cosmology 
with Q m = 0.3, Da = 0.7 and Hq = 70 km s -1 Mpc -1 . 

2. DATA 

The WFC3 Infrared Spectroscopic Parallel Survey 
(WISP) is discussed in full detail in Atek et al. (2010). 
Briefly, WISP consists of HST WFC3 pure-parallel IR 
slitless spectroscopic observations and imaging of hun- 
dreds of uncorrelated high-latitude fields. The spec- 
troscopy is performed using the G102 (0.8 — 1.15 /im, 
R ~ 210) and G141 (1.15 — 1.65 /im, R ^ 130) grisms, 
while the associated near-IR imaging is obtained with 


the FI 10W and F160W, filters. For this paper, we use 
data from 52 separate fields for which both G102 and 
G141 grism spectroscopy are available, covering a total 
of 182 arcmin 2 . These fields include 23 new WISP fields 
in addition to the ones used by Colbert et al. (2013). 
All data are processed with a combination of the WFC3 
pipeline CALWF3 and custom scripts, to account for the 
lack of dithering of the pure parallel data (see Atek et al. 
2010). The siltless extraction package aXe 2.0 (Kiimmel 
et al. 2009) is used to perform the spectral extraction. 
We perform a blind search for emission lines in all fields 
(both grisms) down to a typical 5cr line flux-limit of 
(3 — 5) x 10 -17 erg s -1 cm -2 as explained in Ross et 
al. (2015, in prep). In order to remove the high con- 
tamination rate from false and/or spurious sources due 
to the parallel, slitless nature of the WISP survey, every 
candidate emission line undergoes independent visual in- 
spection by two team members. This process of visually 
confirming the emission lines is described in further de- 
tail in Colbert et al. (2013). 

In this work, we are interested in the Ha A6563 line and 
[Olll] AA4959+5007 doublet, which are covered by WISP 
survey over the redshift ranges: 0.3 < z < 0.7 (only Ha), 
0.7 < z < 1.5 (both Ha+[OlIl]), 1.5 < z < 2.3 (only 
[Olll]). We exclude all sources with any ambiguity in 
their redshift determination among multiple reviewers. 
Specifically, we only retain sources with a quality flag 
< 16, which implies consensus among the independent 
reviewers (Ross et al. 2015, in prep). Since Ha line and 
[Nil] AA6548+6584 doublet is not resolved in the WISP 
grisms, we apply a correction factor of 0.71 to the Ha 
luminosities to account for [Nil] contamination, similar 
to Colbert et al. (2013). Although Villar et al. (2008) 
and Cowie et al. (2011) report decreasing [Nil] /Ha ratio 
with increasing Ha equivalent width, this ratio is nearly 
constant up to Ha EW ~200A, above which the correla- 
tion steepens. The fraction of galaxies in our sample for 
which we may overestimate the [Nil] contribution due to 
the assumed constant correction is only 10%. 

In what follows, we use the completeness analysis from 
Colbert et al. (2013), who performed extensive simula- 
tions to quantify the survey incompleteness as a function 
of line signal-to-noise ratio, equivalent width (EW), as 
well as galaxy size. The completeness simulations fol- 
lowed the full line extraction process after adding arti- 
ficial sources to the real data, spanning a range in red- 
shifts, radii, brightnesses, equivalent widths (EWs), as 
well as using different empirical spectral templates from 
the Kinney- Calzetti Alt as. These simulations show that 
slitless spectroscopic surveys display some level of incom- 
pleteness even at large EWs and line fluxes, primarily 
because of spectral overlap and line mis-identification. 

3. Ha-[OIII] TREND 

Figure 1 shows the observed Ha against [Olll] line lu- 
minosity, for the WISP galaxies in the redshift range 
0.8 < z < 1.2. The red points in Figure 1 show line 
luminosities for a sample of 2141 star- forming galaxies 
in the Sloan Digital Sky Survey (Thomas et al. 2013), 
limited to the redshift range 0.2 < z < 0.3, in order to 
ensure that 3" spectroscopic aperture contains most of 
the galaxy flux (rather than just the central nucleus), 
and to avoid evolutionary effects. 

Despite the large scatter (on the order of 0.5 dex), the 
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Fig. 1. — A broad correlation is observed between H a and [OIII] 
luminosity in the SDSS star-forming galaxies in 0.2 < z < 0.3 ( red 
points ) as well as our WISP sample in redshift range 0.8 < z < 1.2 
( black points). The black and red lines are linear fits to the WISP 
and SDSS data, respectively and the shaded regions show the la 
deviations from the linear fit (~ 0.5 dex). 

two line luminosities are broadly correlated, both in the 
SDSS as well as WISP samples. This is not surprising: 
both Ha and [OIII] are observed in the ionized gas in 
star-forming galaxies, although, while the Ha luminos- 
ity scales directly with the ionizing fluxes of embedded 
young, hot stars, the [OIII] luminosity is more strongly 
dependent on variations in the oxygen excitation state, 
overall gas oxygen abundance, gas density, as well as dust 
reddening (e.g., Kennicutt 1992; Moustakas & Kennicutt 
2006). Figure 1 also shows that the Ha-[OlIl] relation 
does not evolve significantly in the ^4.5 billion years 
elapsed between z ~ 0.25 (SDSS galaxies) and z ~ 1 
(WISP galaxies). Although the gas oxygen abundance is 
observed to evolve with cosmic time, the similarity be- 
tween the observed trends suggests that any evolution is 
masked by the large scatter introduced by the range of 
physical conditions present in galaxies. 

As we will show in Section 7, the observed broad corre- 
lation between Ha and [OIII] luminosity is sufficient for 
the goal of estimating the number of Ha emitters from 
the numbers of [OIII] emitters, as long as the scatter 
is appropriately taken into account. In computing the 
number of Ha emitters from [OIII] emitters at higher 
redshifts, we will assume that the trend between Ha and 
[OIII] does not change in the ^2.2 billion years between 
z ~ 1 and z ~ 1.8. 


(LLF), which describes the volume density of sources as 
a function of both the Ha and [OIII] luminosities. 

The most widely used parametric form for galaxy lu- 
minosity functions is the Schechter function (Schechter 
1976), which is fully described by the parameters L * 
(characteristic luminosity), 0* (number of galaxies per 
unit volume at L*), and a (faint end slope). This func- 
tion is found to reproduce the LLF of both [OIII] and 
Ha-selected WISP galaxies (Colbert et al. 2013). 

We define the bivariate LLF by combining a Schechter 
form (to describe the [OIII] LLF) with the conditional 
probability for finding an Ha source given an [OIII] lu- 
minosity 9 . Thus, we parametrize the [OIII] LLF as: 


^(Loiii) dLom = </>* 



a 

Loin 

{ I* , 

) exp 

L * 


dLom 

L* 


(1) 

where Loin is the [OIII] line luminosity. We adopt a 
log-normal distribution to describe the conditional prob- 
ability that a galaxy with [OIII] luminosity in the range 
(Loin, Lom+dLom) has Ha luminosity in the range 
(LHcn LHa+dLna)- 


?(LHa|L[OIII]) dLua — 

1 

&\nL Ha V2n 


• exp 


ln 2 (L Ha/ (Lhq)) 


2cr lnL H o 


^Lh a 

Lh a 

(2) 


where (Lhq) defines the mean expected Ha luminosity 
for a given [OIII], and cri n L Ha is the scatter around the 
mean relation. (Lhq) and Loin are related through the 
ratio r, such that: 


(Lh a) _ ( Lom^ 

L 0 “ r A Lo ) 


( 3 ) 


The ratio r is defined as the expected Lhcx/Lo at a nom- 
inal luminosity L 0 , where we arbitrarily choose L 0 = 10 4 ° 
ergs s -1 . The LLF and conditional probability equations 
(Equations 1 and 2, respectively) can now be combined 
into the bivariate luminosity function, expressed in terms 
of the Ha and [OIII] logio luminosities, (x and y respec- 
tively) as: 


T(x, y\ P) dx dy = lnlO • 


fio y ' 

a+l 

1(F' 

V 

) exp 



lnlO 

ChnLHc* V^ 


• exp 


2(cr lnLHa /lnlO) 2 


dx dy 


( 4 ) 


4. PARAMETRIZATION OF THE BIVARIATE LINE 
LUMINOSITY FUNCTION 

The goal of this work is to predict the Ha number 
counts at redshift z ~ 2 using the WISP dataset. At this 
redshift, the WISP survey does not cover Ha directly, but 
it does cover the [OIII] AA4959+5007 doublet. Hence, we 
estimate the Ha number counts from the available [OIII] 
number counts. In order to do this, we start by com- 
puting the Ha- [OIII] bivariate line luminosity function 


where (x) is defined by Equation 3 as, 

(x) - log L 0 = log r + fi(y - log L 0 ) (5) 

The bivariate LLF, T(x,^/;P), in Equation 4 can 
now be fully described by the set of parameters P = 

9 In this definition, the marginalized function over all [OIII] 
luminosities is not an exact Schechter form - but very close to it. 
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[a, L*, /3, r, ai n L Hc J. The formulation for the bivariate 
luminosity function described here, has been partly in- 
spired from the size-luminosity bivariate distribution 
from Huang et al. (2013). 

5. FITTING PROCEDURE 

There are various parametric as well as non-parametric 
techniques used to derive the best fit parameters of lu- 
minosity functions (LF); to name a few: the esti- 

mator by Trumpler & Weaver (1953), the C~ method 
by Lynden-Bell (1971), the maximum likelihood estima- 
tor by Sandage et al. (1979, hereafter STY), and the 
stepwise maximum likelihood estimator by Efstathiou et 
al. (1988). In this paper, we use the STY parametric 
maximum likelihood estimator (MLE), modified to ac- 
count for uncertainties in the measurements of the line 
luminosity, as explained in Section 5.2. One of the ma- 
jor advantages of the MLE is that it allows us to fit 
the data without binning. Particularly for small sam- 
ples, this technique reduces the biases introduced by the 
choice of bin-size or bin-center as well as any effects due 
to changing completeness and effective volume within the 
bin (e.g., Mafz Apellaniz & Ubeda 2005). The modifica- 
tion of the method we introduce in Section 5.2 allows us 
to account for significant measurement uncertainties on 
the data. Large photometric uncertainties can impact 
the determination of the best fit parameters of the LLF, 
particularly at the bright end, where the number den- 
sity of galaxies is a steep function of galaxy’s luminosity 
(Henry et al. 2012). In slitless spectroscopy, photomet- 
ric uncertainties can be large even for bright galaxies 
(i.e., the noise is not only due to the sky background but 
also to the possible contamination of the line flux due 
to the continuum of overlapping spectra), so it is crucial 
to account for line luminosity uncertainty in the fitting 
process. 


5.1. Original MLE 

The original MLE is a parametric estimator, where the 
best fit parameters are obtained by maximizing the like- 
lihood function (£) of observing the galaxy sample with 
respect to the parameters of the model. For a given LF 
parametric description T(L), the probability for detect- 
ing a given galaxy with log luminosity L is given by: 


P(Li) 


*(£0 • Vff {Li) 

oo 

J *(L) ■ V eS (L) • dL 

Lum 


(6) 


where V e s(L) is the effective volume of the survey. The 
effective volume varies with the galaxy’s line luminosity 
and redshift, and can be written as: 


Zma x 

V eS (L i ) = / lt^ {z) - C{Li ' z)n{z) ' dz (?) 

Zmin 

where [z m i n , z max ] are the redshift range of the survey, 
dVcomov / dz • dLt is the differential co-moving volume at 
redshift z, C(L, z) is the completeness function, and Q(z) 
is the solid angle covered by the survey. The likelihood 
function for the full sample can then be computed as the 


product of the individual probabilities for all galaxies in 
the sample: 


N 

C = l[P(L i ) (8) 

2=1 

The best fit parameters of the LF can be found by 
maximizing the likelihood function with respect to the 
model parameters. It is mathematically and computa- 
tionally simpler to maximize the log-likelihood function: 


N 

ln£ = yUnP(Lj) (9) 

2=1 

Because this method involves ratios between the differ- 
ential and integrated luminosity functions, the normal- 
ization (</>*) cancels out and, hence, it cannot be deter- 
mined by this likelihood maximization procedure, q 6* can 
be computed following (e.g., Alavi et al. 2014): 

N 

r = (io) 

I V(L)-V eS (L)-dL 

L'lim 

where N is the total number of sources in the sample 
and the survey incompleteness is accounted for by the 
effective volume. 


5.2. Modified MLE 

All astronomical observations have an associated mea- 
surement uncertainty. It is crucial to account for these 
uncertainties, particularly when fitting models that vary 
steeply as a function of the independent variable (e.g., 
at the bright end of the Schechter function). In such 
cases, the best fit-parameters can change significantly, if 
even a few sources are scattered toward or away from the 
bright end due to photometric uncertainties. This prob- 
lem is particularly important for slitless spectroscopic 
data, where line flux uncertainties can be substantial 
even for bright line fluxes, due to the common overlap- 
ping of spectral traces. We modify the original prescrip- 
tion to account for observational uncertainties as follows. 

Instead of calculating the probability that a galaxy 
is exactly at a given luminosity, we marginalize the 
Schechter function over the luminosity error probability 
distribution function, assumed to have a Gaussian form, 
centred at L i: and with standard deviation given by the 
measurement uncertainty cq. In other words, the proba- 
bility P(Li ) that an object has a luminosity L^, given the 
Schechter model, is evaluated by integrating with respect 
to L the convolution of the luminosity function with the 
Gaussian function N{L\{L^ cq}): 
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J ^{L^-V^L^-N^Lua^dL 

P(Li) = 1 ^ s 

J 9(L) • V eS (L) • dL 

Llim 


with, 


N(L\{Li,ai}) 



exp 


({L-Li?\ 1 

V 2of )\ 


( 11 ) 


where L is the log luminosity, V e s(L) is the effective vol- 
ume from Equation 7 - which also accounts for the com- 
pleteness and area coverage of the survey. In the limit 
where the uncertainties are very small, the Gaussian be- 
comes a delta function and the probability approaches 
the value defined in the original MLE, thus, recovering 
the original expression. 

In order to test the performance of our modified MLE, 
we performed a set of simulations (described in Ap- 
pendix A) to reproduce single-line luminosity functions. 
When the sample includes a small number of bright 
sources with significant measurement uncertainties the 
original MLE is less robust than the modified MLE, 
which marginalizes the probabilities over the measure- 
ment uncertainties. For a more detailed discussion, see 
Appendix A. 


5.3. Setting up the Bivariate LLF 

The modified MLE method described in the previous 
Sections can be extended to the bivariate LLF by replac- 
ing the single line LF with the bivariate LLF from Equa- 
tion 4 and marginalizing over both the Ha and [Olll] 
luminosities. The probability for a galaxy with Ha log- 
luminosity (x, in the equations below) and [Olll] log- 
luminosity (y, in the equations below) can then be writ- 
ten as: 


P{xi,Vi) = 


/ 

Llim G) 


y;P) • N{x, y\{xi, {yi, p y ,i})- 


.^^-(z)-C(x i ,z i )-n-dxdy dz 

oo 

J v(x, y; P) • ( z ) ■ c ( x ’ z) -0,-dx dy dz 

Llim G) 

with, 


N(x,y\{xi, a xA }, {yi, a y ,i}) = 

(x-Xif {y-Vi) 2 \ 

Ki Ki J 

(12) 


-exp 


27T (J rp 


where a x and a y are the measurement uncertainties in 
Ha and [Olll] log luminosities, T is now the bivariate 
luminosity function from Equation 4, and C is the com- 
pleteness function (further described in Section 5.4). 


The probability for each source is calculated according 
to Equation 12. We construct the log likelihood func- 
tion as in Equation 9. The log likelihood function is 
maximized and the best fit parameters are obtained us- 
ing scipy . optimize . fmin_l_bfgs_b . Here, all 5 free 

model parameters P = [a, L*, /?, r, <Ji n L Ha ] that define the 
bivariate luminosity function are left free and determined 
by the maximizing likelihood function. 

The scipy . optimize . fmin_l_bf gs_b is a SciPy pack- 
age that uses a limited-memory Broyden-Fletcher- 
Goldfarb-Shanno (BFGS) algorithm in order to find the 
minimum of a function within the parameter space. 
The BFGS algorithm approximates the iterative New- 
ton’s method for finding solutions to functions. The L- 
BFGS algorithm modifies BFGS algorithm to reduce the 
amount of computer memory used and is well suited for 
optimizing functions with large number of variables. The 
version of the algorithm implemented here, L-BFGS-B 
(Zhu et al. 1997), was written by Ciyou Zhu, Richard 
Byrd, and Jorge Nocedal 10 . 

Once the best fit parameters, [a, L*, /?, r, cri n L Hc J 5 are 
obtained, the normalization factor </>* can be computed 
as: 


J y; P) • ( z ) • C{x, z) ■ ndx dy dz 

Llim G) 

(13) 

where the integration limit is taken to be the median flux 
limit for all fields, N is the number of sources detected in 
the survey, x and y are the Ha and [Olll] log luminosi- 
ties respectively, T is the bivariate luminosity function, 
Q is the solid angle surveyed, dV cornov /dz dtl(z) is the 
differential comoving volume at redshift z, and C is the 
completeness function. 

We obtain accurate errors for our best-fit model pa- 
rameters by performing a Markov Chain Monte Carlo 
(MCMC) analysis using the publicly available emcee 
Python package 11 (Foreman-Mackey et al. 2013). We 
use uninformative uniform priors for our parameters and 
the likelihood function is defined as Equation 9 with in- 
dividual probabilities described by Equation 12. 


5.4. Survey Incompleteness 

Colbert et al. (2013) performed an extensive complete- 
ness simulation to quantify the survey incompleteness 
for WISP, which we adopt here. Their completeness is 
provided as a function of equivalent width and signal-to- 
noise ratio. In order to implement the completeness func- 
tion into our formulation, we have to re-par ametrize it 
as a function of flux. The completeness for our sample is 
given by the Ha luminosity (converted to signal-to-noise 
ratio using the survey limit) and marginalized over the 
equivalent width distribution for WISP sources, assum- 
ing that the equivalent width distribution is independent 
of the flux. 


10 http : / / www . ece . northwestern . edu/ ~nocedal/lbf gsb . html 

11 http://dan.iel.fm/emcee/ 
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TABLE 1 

Best Fit Parameters for the Ha-[OIII] bivariate 

LUMINOSITY FUNCTION FIT FOR 0.8 < Z < 1.2 SAMPLE. 


Parameter 

Best-Fit Value 

Oi 

i c +0.5 
-0.2 

log L* 

42 1 +0-1 
-0.2 

log 0* 

-2 05 + 0 - 33 
z.yo _o.l8 

P 

112 +0-06 
- L -- LO -0.26 

r 

0.28 +J5 

ain L Ha 

0.92 ±0 °® 


a CTlnHa = ^logio Ha X ln(10), where C7iog 10 Ha is in dex 


c(i , 2 ) = /c( EW.S/N.^-A^) d(EW) 

(14) 

where L is the line luminosity, EW and S/N are the 
equivalent width and signal-to- noise of the line, fu m is 
the flux limit, and cIl(z ) is the luminosity distance at 
redshift z. 

6. FITTING THE BIVARIATE LLF AT z ~ 1 

In this section, we use WISP galaxies to derive the best 
fit parameters of the bivariate LLF at z ~ 1, where both 
Ha and [Olll] emission lines can be detected in the wave- 
length range covered by the G102+G141 spectroscopy. 
We select the sample to include only galaxies in the red- 
shift range 0.8 < z < 1.2 to allow for sufficient sample 
size while minimizing the impact of an evolving lumi- 
nosity function over the redshift range. We select 487 
galaxies from the 52 WISP fields that satisfy the quality 
flag cut described in Section 2 and have Ha signal-to- 
noise (S/N) ratio > 5. Of this sample, 166 galaxies have 
detected [Olll] with S/N > 2. After applying a strict 
cut of >5 <j in Ha, looking for the [Olll] line is no longer 
a blind search. This allows us to relax the S/N cut for 
the [OKI] line while maintaining a high quality, pristine 
sample. Moreover, we do properly account for the errors 
in the line luminosity during the fitting procedure (see 
Section 5.2). 

We fit the Ha-[OlIl] bivariate LLF to the sample of 166 
galaxies, following the procedure described in Section 5.3. 
The best-fit model parameters are reported in Table 1 
and shown in Figure 2. We also perform the MCMC 
analysis for the model parameters and their posterior 
distributions are shown in Figure 3. 

In Figure 2, the black data points show the sample 
used to fit the bivariate LLF and the contours show the 
best-fit bivariate LLF. The density map shows the Kernel 
Density Estimate (KDE) of the data points (detected in 
both Ha and [Olll] ) corrected for the survey incomplete- 
ness. Both the density map and the contours are plotted 
on the same color-scale. The grey shaded regions show 
the survey flux limits at z = 0.8 (darker) and z = 1.2 
(lighter). 

With the best fit parameters for the bivariate LLF, 
we can derive the single-line luminosity function by 
marginalizing over the nuisance dimension (e.g., the 
[Olll] luminosity function can be computed by integrat- 
ing over the Ha nuisance dimension). In Figure 4, 
we compare the marginalized [Olll] luminosity function 


with the result of the LF of Colbert et al. (2013), com- 
puted for galaxies in the 0.7 < z < 1.5 redshift range. 
Our new estimate of the marginalized [Olll] LF has a 
faint end slope consistent within the errors with the 
slope derived by Colbert et al. (2013, — L5lo'J versus 
— 1.4 zb 0.15). The characteristic T* OIII ] luminosity is 
somewhat lower, although the results are within 2cr from 
each other (42.1^0*2 versus 42.34 ± 0.06). However, we 
note that the two analyses are not expected to provide 
the same best-fit parameters for various reasons: first the 
narrower redshift range used in our work minimizes the 
effect of the evolution of L* over the redshift interval, and 
second our fitting method is not affected by the arbitrary 
choice of bin size and centers, as well as it accounts for 
uncertainties in the line luminosities. 

Now, we can compute the Ha LF by integrating the 
bivariate LLF over the [Olll] dimension and accounting 
for the [Olll] non-detection rate. For fitting the bivari- 
ate LLF, we only used sources detected in both Ha and 
[Olll]. However, there is a significant fraction of sources 
that are detected in Ha but not in [Olll] despite be- 
ing within the wavelength coverage, due to the signif- 
icant intrinsic scatter in the Ha-[OlIl] relation as well 
as the sensitivity limits of our survey. Thus, the bi- 
variate LLF parameters obtained above reproduce the 
number density of galaxies detected in both lines, but 
underestimates the number density of Ha-emitters se- 
lected regardless of their [Olll] luminosity. With the 
goal of obtaining Ha number counts from the [Olll] lumi- 
nosity function, we compute a statistical correction that 
accounts for the fraction of Ha-emitters missed due to 
[Olll] non-detections as a function of the Ha luminosity. 
The non-detection correction term is then applied when 
collapsing the bivariate LLF to obtain the single-line Ha 
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42.5 
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Fig. 2. — WISP data sample for 0.8 < z < 1.2 plotted along with 
the bivariate LLF best fit shown by the contours. The completeness 
corrected Kernel Density Estimate (KDE) map is plotted in color. 
The grey shaded regions represent our survey limit at z = 0.8 
(darker) and z = 1.2 (lighter). 
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Fig. 3. — The posterior as well as joint-posterior distributions 
for the bivariate LLF model parameters obtained from MCMC 
analysis. The best-fit parameters as obtained by the MLE are 
shown by black dots , and for the joint-posterior distributions, the 
68% (thicker) and 95% (thinner) confidence contours are shown. 

LF. 

We compute the non-detection correction term as the 
fraction of galaxies below the detection limit in [Olll] , in 
bins of Ha luminosity. For this analysis we used all galax- 
ies in the redshift range 0.7 < z < 1.5, where both emis- 
sion lines are covered in the spectroscopic observations. 
The non-detection correction ranges between 100% and 
^400%, for Ha luminosities between 10 43 and 10 41 ergs 
s _1 , respectively. In Figure 5, we show the marginal- 
ized Ha LF without the non-detection correction (black 
solid line) and with the correction applied (red line), and 
compare the results with the Ha LF derived in Colbert 
et al. (2013). The marginalized Ha LF corrected for the 
[OIII] non-detection fraction is excellent agreement with 
the Colbert et al. (2013) Ha LF. We note that in our 
formulation, the Ha LF marginalized function over the 
[OIII] dimension is not exactly a Schechter function, but 



Fig. 4. — The collapsed [OIII] LLF derived from the best-fit 
bivariate LLF from Figure 2 compared with the [OIII] luminosity 
function from Colbert et al. (2013). The shaded regions represent 
the la deviations in the best-fit parameters. 



Fig. 5. — The collapsed Ha LLF derived from the best-fit bivari- 
ate LLF from Figure 2 compared with the Ha luminosity function 
from Colbert et al. (2013). The uncorrected (in black) and non- 
detection corrected (in red) LFs are shown. The shaded regions 
represent the la deviations in the best-fit parameters. 

as Figure 5 shows, it is very close to it. 

7. ESTIMATING Ha FROM [OIII] AT z ~ 1 

The main goal of this work is to estimate the Ha num- 
ber counts at z ~ 2, starting from the [OIII] LF at the 
same redshift. Here, we test how accurately the Ha num- 
ber counts obtained with this approximation reproduce 
the known Ha number counts at z ~ 1 , obtained through 
direct integration of the Ha LF. To this aim we use the 
sample of 129 [OIII] -selected galaxies in the 0.8 < z < 1.2 
redshift range and with [OIII] S/N > 5, along with the 
same quality flag cuts as described in Section 2. The 
[OIII] single line luminosity function is computed using 
the modified MLE and the completeness analysis from 
Colbert et al. (2013). The best-fit parameters for the 
[OIII] LF at z ~ 1 are reported in Table 2. 

Figure 6 shows our best-fit [OIII] single line luminos- 
ity function as well as the results from the literature. 
As noted before, we find a good agreement between our 
and Colbert et al. (2013) measurements of the LF on the 
WISP datasets. The rise in the faint-end for Khostovan 
et al. (2015) and Sobral et al. (2015) H/?+[OlIl] LFs can 
be attributed to the Hf3 emitters in their sample. How- 
ever, we note that the variation among different measure- 
ments is still substantial, especially at the bright end - 
the number density of L* galaxies (i.e., ^ 10 42 ergs s -1 ) 
varies by almost an order of magnitude. This comparison 
clearly shows that the nominal errors typically quoted on 
the best-fit LF parameters generally do not provide an 
adequate measurement of the actual variation observed 
in LF determination. The source of this variation is 
most likely systematic (e.g., different selection techniques 
provide systematically brighter /fainter samples, galaxy 
clustering could systematically enhance/suppress num- 
ber counts in small-area fields, and so on), and need to 
be accounted for in predictions used to optimize large 
galaxy redshift surveys aiming to constrain dark energy. 

Using the best-fit Schechter parameters obtained for 
the [OIII] single line LF, we reconstruct the bivariate 
LLF using the parameters /?, r, and cri n Ha from Table 1. 
As done in Section 6, we compute the Ha luminosity 
function by marginalizing over the [OIII] dimension. In 
Figure 7, we compare the Ha single-line luminosity func- 
tion obtained from the [OIII] LF with the direct esti- 
mate from Colbert et al. (2013). Clearly, the two LF 


TABLE 2 

Best Fit Parameters for the [OIII] exclusive 
LUMINOSITY FUNCTION FIT FOR 0.8 < Z < 1.2 SAMPLE. 


Parameter 

Best-Fit Value 

a 

-1 42 + 0 - 23 

- 0.43 

log L* 

42 21 + 0 - 22 
uz.zi _ 0 18 

log 0 * 

q 17 +0.27 
-o.li _ 0 .39 


agree very well, as demonstrated also by the cumulative 
number counts shown in Figure 7, where we show the 
Ha number counts obtained both from the [OlIl]-only 
fit and from the bivariate LLF fit from Section 6 (black 
solid and dashed line, respectively). The errors on the 
number counts account for the uncertainties in the best 
fit parameters in addition to the normal Poisson errors. 
Also note that the number counts are corrected for sur- 
vey incompleteness. As evident from the figure, the re- 
covered Ha number counts from the [OlIl]-only fit agree 
extremely well with the bivariate version as well as from 
direct integration of the Ha LF. We add for completeness 
the number counts from Sobral et al. (2012) and Geach 
et al. (2010). The differences between the WISP dataset 
and these two works are discussed in detail in Colbert et 
al. (2013). 


8. ESTIMATING Ho FROM [OIII] AT z ~ 2 

Having demonstrated the feasibility of our procedure 
at z ~ 1, we now apply it to the redshift 2 case. In Sec- 
tion 3, we compared the Ha-[OlIl] correlation for SDSS 
(0.2 < z < 0.3) and WISP (0.8 < z < 1.2) data. There is 
little evidence for significant evolution of the Ha-[OlIl] 
correlation between z ^ 0.25 and z ~ 1. Continuing 
with the assumption that the Ha-[OlIl] correlation from 
z ~ 1 also holds at z ~ 2, we use the sample of WISP 
[OlIl]-emitters at z ~ 2 together with the Ha-[OlIl] ra- 
tio parameters obtained at z ^ 1 (see Section 6) to derive 
the z ~ 2 Ha number counts. We follow the same steps 
as in Section 7. Namely, we first fit a Schechter model 
to the [OlIl]-only line LF. Next, we use the best-fit pa- 
rameters together with r, <JLi nHa i and (3 from Table 1 to 
construct the z ~ 2 bivariate LLF. Finally, we compute 
the marginalized z ~ 2 Ha LF, and integrate it to obtain 



Fig. 6. — The best-fit [OIII] LLF derived for 0.8 < z < 1.2 
WISP dataset compared with the [OIII] luminosity function from 
Colbert et al. (2013) as well as other estimates from the literature. 
The shaded regions represent the lcr deviations in the best-fit pa- 
rameters. The inset shows the joint-posterior distribution of a and 
L* from MCMC analysis for our best-fit [OIII] LF. 
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Fig. 7. — Top: The collapsed Ha LLF for 0.8 < z < 1.2 derived 
using the best fit [OIII] LLF in Figure 6 and the best fit bivariate 
LLF parameters /3, r, and crinHo; from Table 1. Bottom: The Ha 
number counts estimated for z ~ 1 from our Ha LF, along with 
estimates from other groups in the literature. The errors on our 
estimate accounts for the uncertainties in the best-fit parameters 
in addition to the normal Poisson errors. All number counts are 
corrected for survey incompleteness. 


the Ha number counts. 

The z ~ 2 sample consists of 91 WISP [OIII] -emitting 
galaxies selected to be in the redshift range 1.85 < z < 
2.2, to have [OIII] S/N > 5, and redshift quality flags 
< 16. The [OIII] single line luminosity function is fit us- 
ing the modified MLE, accounting for the measurement 
uncertainties, and the completeness analysis from Col- 
bert et al. (2013). An additional completeness factor is 
applied in order to account for the loss of high - z [OIII] 
emission lines, due to the inability to resolve the dou- 
blet, as discussed in Colbert et al. (2013). The best-fit 
parameters for the [OIII] LF at z ~ 2 are reported in 
Table 3. Our best-fit z ~ 2 [OIII] luminosity function is 
plotted in Figure 8, along with the result from Colbert et 
al. (2013) and other estimates from the literature. Our 
best fit LF shows a slightly steeper faint-end slope and a 
higher </>* than what was derived by Colbert et al. (2013) 
on a smaller sample. The errors in our previous work, 
however, were substantial, and the differences are not 
significant. 

Using the best-fit [OIII] single line luminosity function, 
the bivariate LLF is reconstructed using the parameters 
/?, r, and a\ n n a from Table 1 and the Ha luminosity 
function is obtained by marginalizing over the [OIII] di- 
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TABLE 3 

Best Fit Parameters for the [OIII] exclusive 

LUMINOSITY FUNCTION FIT FOR 1.85 < Z <2.2 SAMPLE. 


Parameter 

Best-Fit Value 

a 

-1.57 

log L* 

42.55 t°;“ 

log b* 

-2.69 t° 0 fi 


mension. In Figure 9, our best-fit z ~ 2 Ha single line 
luminosity function obtained from fitting just the [OIII] 
data is shown alongside other estimates from the liter- 
ature. The variation among different determinations is 
large, probably because of systematic uncertainties due 
to different selection techniques, area covered, and proce- 
dures used for the estimates of the Schechter parameters. 
As noted before, these systematic effects are typically not 
accounted for in the errors quoted alongside the best- 
fit estimates of the parameters. Thus, the shaded areas 
shown in Figure 9, are lower limits to the real variation 
of the volume density at each Ha luminosity. 

9. Ha NUMBER COUNTS 

We use the collapsed Ha luminosity function to com- 
pute the 1.85 < z < 2.2 Ha number counts down to the 
a range of limiting flux values expected to be reached 
by future dark energy surveys, and plot the results in 
Figure 9. For comparison, we have also plotted numbers 
from Geach et al. (2010), Sobral et al. (2013), and Lee 
et al. (2012). The numbers for Geach et al. (2010) are 
reduced by a factor of ln(10) to account for an error in 
the published article, resulting from improper conversion 
of T(log L) luminosity functions to the standard T(L) lu- 
minosity functions. In Figure 9, our number counts are 
higher than previous estimates at all flux limits, although 
below 5x 10 -16 ergs s -1 cm -2 , Geach et al. (2010), Lee et 
al. (2012), and our work agree within the error bars. The 
Sobral et al. (2013) counts are still lower than any previ- 
ous as well as our estimates. At brighter fluxes, the vari- 
ation among the various estimates is very large. Num- 
ber counts of bright rare galaxies, however, are strongly 
affected by sample variance. The WISP number counts 
suffer less from this effect, because of the observing strat- 
egy (52 independent fields scattered over the full sky). 



Fig. 8. — The best-fit [OIII] LLF at z ~ 2 derived from WISP 
data plotted along with other estimates from the literature. The 
shaded regions represent the lcr deviations in the best-fit parame- 
ters. The inset shows the joint-posterior distribution of a and L* 
from MCMC analysis for our best-fit [OIII] LF. 
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Fig. 9. — Top: The collapsed Ha LLF at z ~ 2 derived using the 
best fit [OIII] LLF and the best fit bivariate LLF parameters from 
Table 1. The shaded regions represent the lcr deviations in the 
best-fit parameters. Bottom: The Ha number counts estimated 
for z r-j 2 from our Ha LF, along with estimates from other groups 
in the literature. The errors on our estimate accounts for the uncer- 
tainties in the best-fit parameters in addition to the normal Poisson 
errors. All number counts are corrected for survey incompleteness. 

Finally, we also provide number counts for the whole 
redshift range expected to be covered by the upcoming 
dark energy surveys. Figure 10 shows the expected Ha 
number counts as a function of survey flux limit for the 
redshift range 0.7 < z < 2. The redshift range is broken 
into two: 0.7 < z < 1.5, where the best-fit Ha LF from 
Section 6 is used, and 1.5 < z < 2, where the result 
from Section 8 is used to compute the number counts. 
For comparison, we again plot the Colbert et al. (2013) 
(for 0.7 < z < 1.5) and Geach et al. (2010) (for 0.7 < 
z < 2.0) number count estimates. For the redshift range 
0.7 < z < 2, we expect ^3000 galaxies/deg 2 for a flux 
limit of 3 x 10 -16 ergs s -1 cm -2 (the proposed depth of 
Euclid galaxy redshift survey, see Laureijs et al. (2011)) 
and ^20,000 galaxies/deg 2 for a flux limit of ~ 10 -16 
ergs s -1 cm -2 (the baseline depth of WFIRST galaxy 
redshift survey, see Spergel et al. (2015)), when probing 
with Ha. Number counts for various redshift ranges and 
limiting fluxes are summarized in Table 4. These number 
counts have been corrected for survey incompleteness as 
well as for [Nil] contamination as discussed in Section 2. 

The planned spectral resolution for the Euclid mission, 
at the time of writing, is R/^250 (Laureijs et al. 2011), 
which will not be able to resolve Ha+[NII]. Hence, we 
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also provide the numbers counts that are not corrected 
for the [Nil] contamination - these are summarized in 
Table 5. For the redshift range 0.7 < z < 2, we expect 
^5700 galaxies/deg 2 for the Ha+[NII] flux limit of 3 x 
10 -16 ergs s -1 cm -2 . 

10. SUMMARY AND CONCLUSIONS 

Upcoming space based missions will be performing 
galaxy redshift surveys with the aim of understanding 
the physical origin of dark energy. The constraints that 
a given mission will be able to place on the dark en- 
ergy equation of state parameters depend on the surface 
density of the used tracers. Both Euclid and WFIRST- 
AFTA will be using Ha and [Olll] emitters as tracers of 
the galaxy population and will focus on the 0.7 < z < 2 
redshift range. The precise redshift intervals, however, 
are still being tuned to maximize the scientific output of 
these missions. Here, we use the WISP survey to extend 
on our previous work (focused on Ha number counts up 
to z=1.5) and statistically estimate the number counts of 
Ha emission line galaxies in the full 0.7 < z < 2 redshift 
range. 

To this aim, we have measured the bivariate Ha-[OlIl] 
LLF at z ^ 1, and showed how, at this redshift, Ha num- 
ber counts can be accurately predicted from the [Olll]- 
only line LF, if the relationship between the Ha and 
[OIII] luminosities is known. We find that these two lu- 
minosities are broadly correlated, admittedly with a large 
scatter, that is dominated by different oxygen excitation 
states and amount of galaxy dust extinction. The large 
scatter is observed both in the nearby sample (z ~ 0.25, 
from SDSS observations) as well as at z ^ 1. Moreover, 
we find no significant evolution in the best-fit [OlIl]-Ha 
relation in the ^4.5 billion years elapsed between these 
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Fig. 10. — The Ha number counts estimated for 0.7 <z < 2.0- 
relevant redshift range for future surveys to cover Ha. The redshift 
range for each estimate is reported in parenthesis in the legend. 
Our estimate solid black is split into two ranges: 0.7 < z < 1.5 
(dashed black) and 1.5 < z < 2 (dotted black), which use the z ~ 1 
and z ~ 2 Ha LFs, respectively. We also plot two estimates from 
Sobral et al. (2013): solid green line represents the H a number 
counts estimated using the Ha LF derived at z = 1.47 and dashed 
green line represents the sum of number counts estimated over the 
redshift ranges 0.7-1. 2, 1.2-1.85, 1.85-2.0 using Ha LFs derived 
at z=0.84, 1.47, 2.23, respectively. The errors on our estimate 
accounts for the uncertainties in the best-fit parameters in addition 
to the normal Poisson errors. All number counts are corrected for 
survey incompleteness. 


two epochs. We make the working assumption that the 
relation continues not to evolve significantly out to red- 
shift z ^ 2, or, in other words, that any evolution is 
masked by the large scatter observed in the relation. 

To fit the bivariate LLF model to the data, we intro- 
duced a modified Maximum Likelihood Estimator that 
allows us to properly account for the uncertainties in the 
line flux measurement. This modification can change the 
estimate of the best-fit parameters, particularly for mod- 
els that vary steeply over small range of luminosities. 
Our simulations show that the modified MLE improves 
the accuracy of the recovered best-fit parameters - es- 
pecially, when dealing with larger samples, where the 
measurement uncertainties are more significant than the 
uncertainty introduces by small number statistics. 

We combined the direct measurement of the z ~ 1 — 
1.5 Ha LF with the z ~ 2 Ha LF determined from the 
[OIII] LF and the bivariate LLF information to provide 
an estimate of the number of Ha emitters expected to be 
observed down to different line flux limits. Our number 
count estimates in the full 0.7 < z < 2 are approximately 
40% lower than those of Geach et al. (2010) at the bright 
flux limits (i.e., for line fluxes above 3.5 x 10 -16 ergs 
s -1 cm -2 ), confirming, with twice as many fields and 
with the full redshift range, the result of Colbert et al 
(2013) based on the number of Ha emitters up to 2 ^ 1.5. 
However, we note that the variation in the number counts 
obtained from different published works in the literature 
is substantial at these bright flux levels. This is due to 
a combination of effects, including the different sample 
selection techniques, fitting algorithms used to obtain the 
Schechter parameters, as well the different area/depth 
combinations of various surveys. 

The work and results presented in this paper give us 
a better understanding of the expected performance 
from future planned galaxy redshift surveys aiming 
at constraining the properties of dark energy. This 
is a significant step toward reducing the uncertainty 
of figure-of-merit for dark energy for both Euclid and 
WFIRST-AFTA. In order to further optimize these 
planned surveys, more homogenous data are needed. 
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TABLE 4 

Cumulative Ha number counts (after applying [Nil] correction) a 


Ha Line Flux Limit b 

0.8 < 2 < 1.2 

1.85 < z < 2.2 

0.7 <z < 1.5 

1.5 < z < 2.0 

0.7 < 2: < 2.0 

1.0 

55781“* 

Q7Q0+ 111 

OMU - 1465 

9665l®4g 

9813l^ 6 

19478I2530 

2.0 

2015^263 

976+ 91 
y ‘°-539 

3505+“ 

2848+° 2 

6353l“i 

3.0 

100411“ 

4201**3 

1769l|f 

1282l«« 2 

305212“ 

5.0 

369+“ 

136tgl 

6721S 

440l 2 2 3 6 

11131233 

7.5 

151±& 

59+ 14 

OZ -42 

2851“ 

179+ 14 
1 ' y -106 

464l“ 2 

10.0 

76+1 

25 +1 ° 
zo_ 21 

1 47 + 36 

-43 

91+ 10 

b,± -58 

2381*2 


a All counts are per deg 2 . Fluxes have been corrected for survey incompleteness. The errors account for uncertainties in the best-fit 
parameters along with the normal Poisson errors. 
b in 10 -16 ergs s -1 cm -2 . 


TABLE 5 

Cumulative Ha number counts (without applying [Nil] correction^ 


Ha+[NII] Line Flux Limit b 

0.8 < 2 < 1.2 

1.85 < z < 2.2 

0.7 < z < 1.5 

1.5 < z < 2.0 

0.7 < 2 < 2.0 

1.0 

8614+*g 

67361126! 

150011°“ 

169241““ 

319251““ 

2.0 

34171 233 

19221“° 

59171“° 

5365l°298 

11282+116 

3.0 

i8i7iiii 

8571*11 

31661“° 

25241*“ 

5690+“ 

5.0 

700+95 
' °°-258 

9Q1 + 140 

-169 

13031“! 

908li°° 

991 0+243 
^J-Z _ 456 

7.5 

qoq+67 

ozo -93 

H7 + 78 
11 ' -77 

Soling 

384l“* 

Q75+151 
y ‘ °-234 

10.0 

171+ 37 

4 -58 

5Q+52 

oy -43 

3211to 

2021®*, 

504+ 93 

O^T_i4o 


a All counts are per deg 2 . Fluxes have been corrected for survey incompleteness. The errors account for uncertainties in the best-fit 
parameters along with the normal Poisson errors. 
b in 10 -16 ergs s -1 cm -2 . 
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APPENDIX 

A: I D SIMULATIONS TO TEST THE MODIFIED MLE 

Before applying the modified MLE to deriving the H<a-[OlIl] bivariate LLF for our 0.8 < z <1.2 WISP sample, we 
test the validity of our modifications to the MLE - results of which are expected to scale to the bivariate case. 

We generate 1000 samples of galaxies distributed according to a known luminosity function. We run two sets of 
simulations for two different sample sizes: (i) small (200 sources per sample), roughly the number of sources in our 
sample, and (ii) large (2000 sources per sample), roughly the number of sources expected to be covered by the end of 
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the WISP survey within the redshift range of interest. The simulated galaxies are assigned the typical uncertainties 
observed for WISP galaxies at similar luminosities are further randomized within that error-bar. We then fit the 
simulate samples with both the original and modified MLE techniques. 

Figure 11 shows the results for the single line LF simulations, for the two different sample sizes for both the original 
and modified MLE. The modified MLE recovers the true parameters with greater overall accuracy in both large and 
small sample size cases, even though the scatter is similar. Since measurement uncertainties are not properly treated 
by the original MLE, a few bright sources with large uncertainties can skew the results significantly. The modified 
MLE is much less prone to this effect since it marginalizes over the uncertainty. The efficiency of the two estimators 
depends on what factor is dominating: the statistical randomness of the sample or the measurement uncertainties of 
the sample. 

Since WISP is a slitless grism spectroscopy survey, even the bright sources can have significant uncertainties due to 
crowding, contamination or other issues. For our sample, the modified MLE is expected to provide an improvement 
over the original MLE. 



Fig. 11. — Results from the 1-D simulations for small (200 sources per sample; top row) and large (2000 sources per sample; bottom 
row) sample sizes comparing the original ( left column) and modified ( right column) MLEs. The solid black lines show the true parameters 
expected to be recovered by the fitting procedures. The contours are at 10%, 33%, 66%, and 95% confidence levels. 


